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This note describes a random walk equivalent to the Neumann series solution of an 
integral equation. A diffusion analogy and the problem of importance sampling are 
discussed briefly. 



Certain types of matrices can be inverted by a 
Monte Carlo method that was devised by J. von 
Neumann and S. Al. Ulani, and that appears in a 
paper written by G. E. Fors^^the and R. A. Leibler.^ 
In the following it is shown that this metliod, when 
suitably generalized, can be used to approximate the 
Neumann series solution of an integral equation. 
This procedure is then shown to be similar to those 
used to solve cert am problems debating with scattering 
of particles, which can l)e formulated in lei-nis of a 
nonhomogeneous integral equation, and wliicli have 
been described by llei'maii Kahn.'^ 

We wish to solve the following equation for </> (.t) 



(i>{x)=j{x) + \\ K{x,y)ck(y)dy, 



(1) 



where we assume, for convenience, that K{x,,y) is 
uniformly continuous and a and h are finite. A 
formal solution of (1), which is valid if X is smaller 
in absolute value than the eigenvalue of the cor- 
responding homogeneous equation, which is smallest 
in absolute value, is given by Neumann's series ^ 

cl>{x)=J{x) + \j'K(x,y)f{y)dy 

+ \'j'K^H^,y\ny)dy+-", (2) 

where 

K'^Kx,y)= rK^''-'Hx,OK{^,y)d^. 

J a 

"We may write this as 

0(x)=/(x) + xj L(x,y)j'{y)dy, 
where 

L{x,y)=K{x,y) + \K^'Hx, y) + \'K^'\x,y)+ • • -. 

We need consider hereafter only the series (2) — if 
we wish to obtain the resolvent kernel L, we may 
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use the fact that L satisfies the following equation 

L(x,y)=K{x,y) + \ V K{x,QUi,y)dl 

J a 

The series (2) can be evaluated by a Monte Carlo 
procedure consisting of two steps: first, a random 
procc^ss by which a particular term of the series is 
selected, and second, a random process by which the 
multiple integral defining the term selected is ap- 
proximated (see footnote 3). This, when multi])li(Hl 
by a suitable weighting factor, will ])rovide an 
estimate of the vahie of </>(./:) for a ])articular x. We 
notice, howevcM', that in series (2) each term is ob- 
tain(Hl by a sim])le operation on the preceding one, 
and this allows us to consti'uct a simple i-andom 
walk for apj)roximating (2). 

This random walk is defincMl by a sequence of 
functions Pn {x,y) , defined for a ^xkb, a^ y ^b, and 
ri=:0, 1, 2, . . . , where Pn{x,y) >0 and 

£Pn{x,y)dy<l-8, 5>0. 
In terms of this seciuencc* we also deiine 



and 



V (x v) = ^^^\ 
'Pn{x)=l— j Pn{x,y)dy. 



In addition, we must also require that |X| be smaller 
than the smallest eigenvalue of K^{x,y) = \K(x,y)\, 
To estimate 0(x), we start a random walk at the 
point x of the interval {a,b). With a probability 
density given by Po(x,Xi) we jump to the point Xi, or 
we stop the random walk, with probability po(x). 
Suppose that we have made k jumps, and have 
reached the point x^. Then with probability density 
P}c(Xk,x,c^i) we jump to the point j^a:+i, or we stop the 
random walk with probability Pkixjc)- If our ran- 
dom walk has traversed the path y:x—^Xi-^ . . . -^ 
Xjc-i-^Xjc, stopping at Xjc, we define 



V(y/x) = Vo{x,Xi) . . .Vk-i(Xk-i,Xi) 



pjcixy 



(3) 



and use Viy/x) as an estimate of <^(x). The expecta- 
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tion of V(ylx) is defined by 

E{V(y[x))== f V(y/x)dP{ylx)- (4) 

J Ally 

Since the probabilities used at each step are con- 
ditional probabilities, the probability measure de- 
fined in the space of paths with exactly k steps is 
given by 

dP(y/x)=Po(x,x,) . . . 

Pjc-l(Xk-i,Xjc)p!c(x,c)dXi . . . dXfc, (5) 

and thus the integral (4) can be represented b}^ a 
sum of integrals, the kih of which corresponds to 
those 7 with k steps. In this way we obtain 



E(V(y/x))=f{x) + \ rK(x,x,)J(xOdxi 

J a 

nb 
K{x,Xi)K(xi,X2)J(x2)dxidx2 + . 

= (l){x), 



(6) 



which proves that V(y/x) is an unbiased estimate of 
(pix). Because of the condition on y, the above in- 
tegrals and series are absolutely convergent. 

Similarly we obtain for the nth moment of the 
V{y/x), (assuming the nth moment exists), 



..(x)=£;(F%/x))=^3^ + Xj^ K(x,xOFS-^^ 

+ X2 K {x, Xi)K {Xi, X2)Vr^ {x, XiW'i-^ {xi, X2) -in-ui\ ^^2 dxi 

J a r2 V^2J 



(7) 



In the particular case that Pn(x,y) is independent 
of n, that is the probabilities depend only on where 
we are, and not on how many steps we have taken, 
(7) reduces to a rather simple form. In fact, if we 
define Kn{x,y)=K{x,y)V''~^{x,y) the nth moment 
satisfies the following equation 



Vn{x) = 



JHx) 



p^'-^ix) 



^X 



X Kn{x,y)vn{y) dy, (8) 



(assuming | X | is sufficiently small) . 

Let us now turn our attention briefly to the prob- 
lem of finding a set of functions Pn(x,y) such that 
the variance of the estimate of (l>(x) will be as small 
as possible. We might expect that if we knew (l)(x) 
we could find a Monte Carlo procedure that would 
give us the correct answer with probability 1 (see 
footnote 3). It is convenient to use an argument 
based on the two-step method mentioned earlier. 
In the series (2) denote the nth term by In{x). 
Our Monte Carlo process gives us a set of functions 
P(n,x) and v(n,x); where P{n,x)v{n,x)=In{x), and 
P{n,x) is the probability of taking exactly n steps if 
we start at the point x. Referring to (5), we see that 



P{n,x)-- 

b 



J a 



Po{X,X\) ' ' • r n-\\Xn-i,Xn) Pn\X n) dXy • • • dXri' 



We would like to choose P{n,x) so that P{n,x)=: 
In(:x)l(t>(x), as then v{n^x) = ^{x). To estimate Inix) 



correctly with probability 1, we must have K(x,y) 
and f{x) non-negative, as we cannot allow negative 
probabilities. It may be convenient to introduce 
a sequence [Unix)}, which has a known sum S'(x), 
and estimate <l)(x)-\-S{x) by the quantity v*(n,x) = 
v(n,x) + Un{x)IP{n,x). 

The above discussion is of course of little more 
than academic interest, but it gives an insight into 
the problems involved in minimizing the variance 
of the estimate. One thing we notice in particular 
is that the optimum choice of the sequence {Pn(x,y)} 
will depend, in general, on the point x at which we 
wish to find </>. 

As promised in the introduction, this method 
can be shown equivalent to Monte Carlo techniques 
often used to solve problems formulated in a 
different manner. Suppose we consider the prob- 
lem, important in certain apphcations, of the diffusion 
of the neutrons through a thick slab of scattering 
material. Let the slab lie between the planes x=^0 
and x^a. Let /(£', 6, cj), x) give the number of 
neutrons incident with energy E, at an angle given 
by d and 0, which make their first collision at a 
distance x from the front of the slab. Let 
K{E\ 6\ <^^ x\ E, 6, (j), x) give the probabilit}^ that a 
neutron that makes a collision a distance x from the 
front of the slab, and has an initial energy E and an 
initial direction of travel (^, <^), will proceed in the 
direction {d\ 4>'), have energy k\ and make its next 
collision at a distance x' from the front. A colli- 
sion at x>a corresponds to a neutron passing through 
the slab, and a collision at x<CO corresponds to a 
neutron reflected back, which we assume stops the 
process. Let yp{E^ 6, <^) give the number of neutrons 
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^^ that penetrate the slab with energy E, and have 
the direction of motion 6^ <t>. 
Then clearly 



^ 



^l^{E,^,<i>)= r f{E,e,<j>,x)dx 

J a 

+ K{E,d,4>,x,E',B',<i>',x') 

Jx+a Jx'=0 jE',e',4>' 

' J(E', e\ ct>', x')dE'dd'd<t>'dx'dx (9) 

+ higher terms representing multiple collisions. 
This is very similar to equatiop (2). 

The usual Monte Carlo method for evaluating 
(9) is to consider a randomly selected incident particle 



and trace its path thru the slab, using some chance 
device to determine the result of a collision (see 
footnote 3). The method we have found for 
evaluating series (2) corresponds to choosing a par- 
ticular particle that has penetrated the slab, tracing 
its path backward thru the slab, and then muUi- 
plying by the density of particles that have the same 
entrance characteristics. Although both the posi- 
tion of the particle making the random walk and its 
final weight must be remembered from step to step, 
only the final weight is needed for the answer in 
our method, while in the other both the final position 
and final weight are needed, and a histogram 
must be plotted. Thus if an answer is desired for 
only a few points, our method may be simpler. 

Los Angeles, August 16, 1950. 
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